Dynamics of a Limit Cycle Oscillator under Time Delayed Linear and Nonlinear 

Feedbacks 
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We study the effects of time delayed Unear and nonhnear feedbacks on the dynamics of a single 
Hopf bifurcation oscillator. Our numerical and analytic investigations reveal a host of complex 
temporal phenomena such as phase slips, frequency suppression, multiple periodic states and chaos. 
Such phenomena are frequently observed in the collective behavior of a large number of coupled 
limit cycle oscillators. Our time delayed feedback model offers a simple paradigm for obtaining 
and investigating these temporal states in a single oscillator. We construct a detailed bifurcation 
diagram of the oscillator as a function of the time delay parameter and the driving strengths of the 
feedback terms. We find some new states in the presence of the quadratic nonlinear feedback term 
with interesting characteristics like birhythmicity, phase reversals, radial trapping, phase jumps, 
and spiraling patterns in the amplitude space. Our results may find useful applications in physical, 
chemical or biological systems. 
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I. INTRODUCTION 

Coupled limit cycle oscillators have been extensively 
studied in recent times as a mathematical model for un- 
derstanding the collective behavior of a wide variety of 
physical, chemical and biological problems p]-p^. One of 
the simplest and earliest of such models is the so called 
Kuramoto model which is a mean field model of a 
collection of phase oscillators, and clearly exhibits such 
cooperative phenomenon as spontaneous synchronization 
of the oscillators beyond a certain coupling strength. A 
more generalized version of the coupled oscillator model 
that includes both phase and amplitude variations ex- 
hibits collective behavior like amplitude death, where, for 
a large enough spread in the natural frequencies of the 
oscillators, an increase in the coupling strength induces a 
stabilization of the origin, leading to a total cessation of 
oscillations in the system . Other collective states 

observed in these models include partial synchronization, 
phase trapping, large amplitude Hopf oscillations and 
even chaotic behavior |]l5|, 16|. Recently there has been 
some interest in investigating the effect of time delay on 
the collective dynamics of these coupled models [p7|-p4|. 
Time delay is ubiquitous in most physical pa- BSl, c hem- 
ical [|9|, biological [Q, neural ecological [p2[, and 
other natural systems due to finite propagation speeds 
of signals, finite processing times in synapses, and finite 
reaction times. Time delayed coupling introduces inter- 
esting new features in the collective dynamics, e.g. simul- 
taneous existence of several different synchronized states 
20[ , regions of amplitude death even among identical 



oscillators |^,g2|, and bistability between synchronized 
and incoherent states p0| , p3[ . 

One of the remarkable aspects of this cooperative dy- 
namics is that many of its salient features can be observed 
even in a system consisting of just two coupled oscilla- 
tors |l|,|l|,|l|,|2|. The temporal behavior of either of the 
two oscillators in such a case (which is easy to investi- 
gate both numerically and analytically) reveals a great 
deal about the collective aspects of larger systems. In 
fact, a useful point of view to adopt is to regard each 
oscillator as being driven autonomously by a source term 
that represents the collective feedback of the rest of the 
system. Motivated by such a qualitative consideration, 
we have studied in detail the dynamics of the following 
model system of an autonomously driven single limit cy- 
cle oscillator 



Z{t) - (a + tuj- I Z{t) \')Z{t) = f{Z{t - r)), (1) 

where Z(t) = X + iY is a complex quantity, to the fre- 
quency of oscillation, a a real constant, and r > is the 
time delay of the autonomous feedback term /. In the 
absence of the feedback term Eq. (]l|), often called the 
Stuart-Landau equation, has a stable limit cycle of am- 
plitude ^/a with angular frequency lu. It is simply the 
normal form of a supercritical Hopf bifurcation and is a 
useful nonlinear model for a variety of physical, chemi- 
cal and biological systems. For the autonomous feedback 
term we choose the following model form: 



f{Z{t - r)) = -KiZ{t - r) - K^Z^t - r). 



(2) 
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where Ki and K2 represent the strengths of the Hnear 
and nonhnear contributions of the feedback. This choice 
is motivated by considerations of both mathematical sim- 
pUcity and possible importance for modeling of physical 
and biological systems. The quadratic term is the sim- 
plest nonlinearity that can break the rotational symme- 
try of the Stuart-Landau system. Physically this term 
introduces nonlinear mode coupling, a process that is 
important in large coupled systems. Eq. can also 
be viewed as a prototype equation arising in the delayed 
feedback control of an individual physical or biological 
entity that can be modeled by the normal form. Our re- 
sults may thus be of more general and direct utility in 
addition to provide useful insights into the collective dy- 
namics of large systems. Similar studies (using a variety 
of feedback terms) exist for the damped harmonic oscilla- 
tor, e. g. |]33|] , but we are not aware of such investigations 
for our model limit cycle oscillator. A few investigations 
in the past have restricted themselves to the study of 
noise and perturbations ||3^-|3^ on the dynamics of such 
an oscillator. 

The organization of our paper is as follows. In Section 
U, we analyze the dynamics of the oscillator using just 
the linear feedback term and discuss the analytic condi- 
tions for the stability of the origin and the existence of 
periodic orbits. Detailed bifurcation diagrams are plot- 
ted as a function of the various system parameters like a, 
K and t and their similarity to collective states of larger 
systems is pointed out. We also present numerical re- 
sults on higher frequency states, which can coexist with 
the lowest periodic state and discuss the phenomenon of 
frequency suppression of these states as a function of the 
time delay parameter. Section [I] treats the full feedback 



term by including the quadratic nonlinear contribution. 
The bifurcation diagram is a great deal richer now due to 
the existence of two other equilibrium points in addition 
to the origin. We analyze the stability of these equilibria 
and the consequent temporal behavior of the oscillator in 
various parametric regimes. Some novel temporal states 
are pointed out. Section IV summarizes our results and 



discusses their significance and possible applications. 



II. TIME DELAYED LINEAR FEEDBACK 

We begin our analysis of the model Eq. (|l|) by consid- 
ering only the linear feedback term (i.e. K2 = 0), so that 
we have 

Z(t) = (a -I- iLu- I Z{t) \^)Z{t) - KZ{t ~ r), (3) 

where we have put Ki — K for simplicity of notation. 
Note that the above linear feedback term is similar in 
form to the feedback term used extensively in experi- 
mental and theoretical investigations of control of chaos 
using the Pyragas method |3^. The actual form in the 
Pyragas method is [Z{t—T) — Z{t)], which is equivalent to 
replacing the constant a by a -I- if in the above equation. 
However, unlike the systems investigated for the Pyragas 



method, Eq. (g) has no regimes of chaotic behavior. We 
will examine instead the effect of the time delayed feed- 
back on the stability of the origin and on the nature of 
the periodic solutions. 

In the absence of time delay, it is clear from inspec- 
tion that Eq. (^ has a time-asymptotic periodic solution 
given by Z{t) = y/a - ife*"* for a> K. If a < then 
the origin is the only stable solution; i.e., no oscillatory 
time-asymptotic solutions are possible. At a ~ K, the 
oscillator undergoes a supercritical Hopf bifurcation. We 
now analyze systematically the effect of time delay on 
the stability of the origin and the periodic solutions. 



A. Stability of the origin 

The origin Zp — (0,0) is a fixed point of Eq. (||). 
To study its stability we assume that the perturbations 
about Zp grow as e^*, where A is a complex number. Sub- 
stituting in Eq. (|^) and linearizing about Z = Zp, one 
easily obtains the following characteristic equation: 



X = a±iu! — Ke 



(4) 



where the ± sign arises from considering the complex 
conjugate of Eq. (||). This ensures that we have the 
complete set of eigenvalues. For r = 0, one obtains 
A — a — KzLi Lu. The origin is stable in the region of para- 
metric space where Re(A) < 0, which occurs when K > a. 
So the critical, or the marginal stability curve is given in 
this case hy K ^ a. When t ^ 0, Eq. (|) remains a 
transcendental equation with a principal term A e^'^ and 
hence the equation possesses an infinite number of com- 
plex solutions. Let these roots be ordered according to 
the magnitude of the real parts: {. . . , A„, . . . , A2, Ai, Aq}, 
where Re(Am-i) < R.e(Am). The problem of finding the 
stability criterion then reduces to that of finding the con- 
ditions on K, oj and t such that Re(Aj) < 0, for all j. Let 
A = a -|- if3, where a and f3 are real. By substituting this 
in Eq. (^, we get 



a~ a — K e "'^ cos(/3t), 
(3=±uj + K e""^ sin(/3T) 



(5) 
(6) 



We can arrive at the following two equations for a and (3 
by squaring and adding the above two equations and by 
dividing the first equation by the second respectively: 



a= a — [fi — w)/ tan(/3r), 



a- 



a)2, 



(7) 
(8) 



where, in the above and from here after, we consider only 
one set of curves by choosing [i ~ +ti> ± ^77 . The other 
set of curves arising due to (3 = —u ± ^77 is implicit in 
the above since the eigenvalues always occur in complex 
conjugate pairs. From Eq. (Q) we see that (3 is real only 
when {a — a)^ e^""^ < . So for any finite value of K, 
the value of a is bounded from above. 
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FIG. 1. The stability of the origin is shown in a (r, A'"^) diagram for (a) lj = l,a = 1 and (b) lo — 30, a = 1. At lower 
values of there is only one death region that is confined between r = and t = t\ (0, if), but for higher values of lo there are 
generally multiply connected regions of amplitude death. The curves a„ and 6,i represent t\ (n, K') and T2 (n, K') respectively. 
The dotted horizontal line is Ij fiio). 



To obtain the critical curves, set a = 0. This gives 



(9) 



By inverting Eq. (^), and noting that sin(/3+r) > 0, 
sin(/3_r) < 0, and that can be either positive or neg- 
ative, we obtain the following two sets of critical curves: 



Ti(n,if) = 



Inn 



-\a/K) 



bj + \/K'^ - a? 
2mT — cos^^ [a / K) 



(10) 
(11) 



Lo — V — 

In Eq. (|l|), n = 0, 1, 2, . . .. In Eq. (|l]), if /3_ > 0, n = 
1,2,...; if /3_ < 0, n = 0,1,2,.... Thus the critical 
curves exist only in the region K > a. Since, for t = 0, 
the region of stability of the origin is given by K > a, the 
corresponding region, for r > 0, will be given by the area 
between r = and the critical curve closest to the line 
r = 0. This critical curve should be the one on which 
da/dr > 0. From (|) 



dX 
d7 



,-At 



1 - Kre 



-At 



(12) 



and 



da 



Re 



K(i/3)e- 



-i/3r 



1 - Ktc-'P-' 



(3 KsmiPr) D' 
f3{P^uj) D-\ 



(13) 



where D = [1 
real. Hence 



da 
d^ 



K T cos{ fir)]'^ + [Kt sm{PT)]'^ is positive 



> on n, 

> on T2ii K > /(w), 
= onra if if = /(w), 
< on T2ii K < f {uj), 



(14) 



where f{uj) — Va^ + w^. The above condition implies 
that there can be only one stability region if if > /(w). 
There is a possibility of multiple stability regions if 
K < f{uj). Our numerical plot in Fig. |^(a) of the curves 
Ti{n,K) and T2{n,K) reveals that the region between 
T — Q and t — ri(0, K) is the only stability region pos- 
sible for small values of oj. However, as the value of 
LO is increased, the stability regions can be specified by 
< T < Ti(0,if) and T2{n,K) < t < Ti{n,K), where 
n > 0. In Fig. Wo) the critical curves are plotted from 
Eqs. (|l^) and ( |1 1| ) for such a large value of w, namely 
uj = 30, and the multiple stability regions are represented 
by the shaded portions. Note that a similar situation 
arises in the case of two or more limit cycle oscillators 
that are coupled by a time delay. This was investigated 
in detail in H^,^, where the collective stability regions 
were termed amplitude death regions or death islands in 
the K — T space. 

We conclude this section by carrying out a stability 
analysis of the origin in the (a, K) plane for a fixed value 
of r. Using Eqs.(|5|-|6|) we write below the critical curves 
which are non intersecting: 



sin(/3T) 



and a{P) 



l3±uo 
sin(/3T) 



cos(/3r). (15) 



We note that the above expressions for if (/3) and a(/3) 
have singularities aX (3 — n^n/r and between any two 
successive singular points the expressions produce con- 
tinuous curves in (a, K) plane. Following Diekmann et 
al. 1^^, we define the following intervals where the sign 
of the superscript of / indicates the sign of the function 
sin(/3r) in that interval: 
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In - ((2n - 1) J, 2n^) , /+ = (2n^, (2n + 1)^) , 

(16) 

for n = 0,l,2,.... We restrict our attention, without loss 
of generality, to the case of /3 > 0. Hence we can define 
the following curves in (a, K) plane. 



[a,K) 



iP-Lu) cos(/3r) (/3-w) 
sin(/3r) ' sin(/3r) 




ia,K) 



il3 + uj) cos(/3r) (/3 + w) 



sin(/3r) 



sin(/3r) 




These curves are parametrized by (3. The curves are de- 
generate at w = WK It. For w = in-K /t, n = 0, 1, 2, . . ., 
both the curves C± and merge, and there is another 
curve in addition to the above, defined by 



CR^{{K,a) \ a^K} 



(19) 



with f3 = 27171 /t. Likewise for lo = (2n + l)7r/r, n = 
0, 1,2,.. ., the corresponding additional curve is defined 

by 



CR^{{K,a) I a^~K} 



(20) 



with (3 = {2n + 1)ti/t. We are now only left with 
the task of finding out the number of the eigenvalues 
in the right half plane on either side of the curves 
and D^. For our particular problem it is possi- 
ble to carry out this analysis in an exact manner. Let 
F(a, K,X) — \ ~ iuo + Ke^'^'^ = be the eigenvalue 
equation. Define Gi =ReF, G2 — Im F, and at a = 
define a matrix M by 



M = 



dGi 
da 

dG2 
da 



dGi 
dK 

dG2 
dK 



(21) 



We now make use of a proposition of Dickmann et al. 
[ p8[ to determine the positions of the eigenvalues in the 
complex plane with respect to the curves {a{(3) , K {(3)) . 
The proposition states that the critical roots are in the 
right half-plane in the parameter region to the left of the 
curve {a{j3) , K {(3)) , when we follow this curve in the di- 
rection of increasing (3, whenever det M < and to the 
right when det Af > 0. In the present case, on both the 
curves and D^, the matrix is given by 
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FIG. 2. Critical curves in (a,K) plane. The numbers indi- 
cate the number of eigenvalues in the right half plane of the 
complex eigenvalue space. The amplitude death region is the 
region where there are zero eigenvalues in the right half plane. 
The arrow on each curve shows the direction along which /3 
increases. 

and hence det(M) = sin(/3T). It is easily seen that 
det M > on C+, 1)+, and det M < on C;;, D'. And 
finally we know that the stability region for t = is given 
hy a < K. Combining these facts we present our results 
for the stability region and the number of eigenvalues in 
Fig. ||. The amplitude death region is the region where 
there are no (i.e. zero) eigenvalues in the right half plane. 



B. Periodic Solutions 



We now examine the region where the origin is unsta- 
ble. For T = this region sustains periodic solutions 
as discussed in the introductory remarks of this section. 
We now look for periodic solutions in the presence of time 
delay. For this, it is convenient to cast Eq. (^) in polar 
form. For simplicity, we also set a = 1 implying that the 
oscillator without any kind of feedback has a unit circle 
as its periodic solution, and the phase increases linearly 
on the circle. Writing Eq. (^) in polar coordinates we 
have. 



fit) = [1 ~ r'^{t)\r{t) - Kr(t - r) cos[6l(t - r) 



mi 

(23) 



M 




(^2) d{t) =uj~k ''^^ sin[g(t - r) - d{t)\. (24) 
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Recall that when r = 0, Eq. 
totic periodic solution Z{t) - 



has the time asymp- 
lujt non 



VT^n? e 

zero T we can still assume a periodic solution of the form 
Z{t) = R e*^*; i.e., we are looking for solutions of the 
form r{t) — R and 9{t) — ilt, where R and fl are real 
constants. In fact, it can be checked that this is the only 
solution which has a linear growth of the phase. Substi- 
tuting this form in Eqs. (|23| ) and (24) and after some 
algebra, we obtain the following relations for the ampli- 
tude and the frequency of the oscillator: 



i?= x/l-K cos(17t), 
n=oj + K sin(ilr). 



(25) 
(26) 



The oscillator can now lie outside the unit circle when 
cos(r2r) < 0; that is, the amplitude of the limit cycle can 
increase beyond unity if r e {{2n+^)TT/n, (2n-|- |)7r/r2). 
However, the amplitude is bounded for any value of 
r because max{R) < Vl + K. So for any given K, 
the amplitude of the limit cycle stays in the interval 
[^1 — K, ^/l + K] for arbitrary values of r. The above 
condition on fir can be used in Eq. (Eq) to infer bounds 



on the frequency of the limit cycle: uj-'K<rt<uj + K. 

Eq. ( p6| ) admits multiple solutions for the frequency fl. 
The left-hand side, yi = fi, of ( [26| ) is a straight line and 
the right-hand side, 1/2 = uj + K sin(f2T), is a sinusoidal 
curve with amplitude K and a shift of u above the hor- 
izontal axis, y = 0- The multiple solutions (frequencies) 
are given by the intersection of the curves yi and y2- As 
the value of K is increased for a fixed value of cj, the curve 
?/2 makes more and more intersections with yi and thus a 
set of multiple frequencies comes into existence. Similarly 
a variation in lu will also bring about changes in the num- 
ber of possible solutions. But for all values of lo, including 
for u! — 0, multiple frequency solutions are possible. The 
existence of multiple frequencies is a characteristic fea- 
ture of time delay systems and has been noted before in 
the context of the Kuramoto model with time delay jl^] 
and in other studies ||l7,^,^. These multiple frequency 
states that coexist with the lowest frequency state can 
be accessed by a suitable choice of initial conditions and 
have potential applications in coupled oscillator models 
of the human brain. 
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FIG. 3. Finite time delay leads to the existence of multiple periodic states, which can coexist with the primary limit cycle 
state. The frequencies and amplitudes of these states (for r = 0.68) are plotted as a function of the coupling strength K in (a) 
and (b) respectively. As the value of the strength of the feedback, K, is increased, the number of periodic orbits increases. In 
both the plots the dashed portions of the curves are unstable regions. 



In Fig. ||(a) we plot these multiple frequencies, 17, 
for LU ^ 30 and / = 3.25 where / = ujt/{2tt). Fig. 
^(b) shows the corresponding amplitudes of the multi- 
ple states. Some of these states merge if / is an integer. 
This can be inferred from the intersections of the ampli- 
tude curves in Fig. ^(b). At these points, the oscillator 
has two frequencies with a single amplitude. To find out 
these frequencies, let fii and be the two frequencies 
at these degenerate points. Substituting these values in 
the expression for i? in (25), we get ^ = ^ -f- y , where 
m is an integer. 



Stability of the Multiple States. The stability of these 
multiple periodic solutions can be obtained by lineariz- 
ing about each of the solutions. The linearized matrix of 
Eqs. ( p3| ) and (|2j) about the periodic solutions can be 
written as 



Ml = 



A- 



RC{1 



-At 



^) B{l-e-^^) 



where A = I - ZR^, B = K cos(rjr), and C = 
K sin(riT). The corresponding eigenvalue equation is 



5 



det(Mi - XI) = 0, or 



Aie 



-2At 



-A 



26 



A3 + - A4X + A^Xe-^^ = 0, 



where Ai ^ B'^ + C^, A2 = AB + - 2C^, A3 = 
AB + C^, A4 = A + B and A5 = 2B. If r = 0, then 
A — {0, —2(1 — K)}. We thus recover the result that the 
periodic solution is stable for r = when K < 1. For 
T ^ we need to solve the eigenvalue equation numeri- 
cally. Our numerical results for the stability of the mul- 
tiple periodic states are incorporated in Fig. || and Fig. ^ 
where the dashed portions of the curves indicate unstable 
regions. At large value of r, the frequency of the oscilla- 
tion , fi, gets reduced. This is true for all the multiple 
frequency states that the system possesses. In Fig. ^(a) 
the normalized frequency of each of the states is plotted 
against the time delay on a log scale. The frequencies are 
suppressed at a rate proportional to r~^. Fig. ^(b) shows 
the corresponding amplitudes plotted against time delay. 
This feature of frequency suppression has been observed 
in the past for large coupled systems and once again 
seems to have its roots in the behavior of a single oscil- 
lator in the presence of a time delayed feedback drive. 
However, for short time delays, as may be seen from the 
figures, the effect of time delay is some what distinct - 
both the frequency and the amplitude of the oscillator 
increase slightly with r. To understand this behavior we 
present in Section II C analytic expressions for the time 



evolution of Z{t) in the limit of small time delay. 



C. Small r approximation 

For short time delays, it is still worthwhile to expand 
the delay variable Z(t — r) in a Taylor series despite oc- 
casional warnings |39| ]. Let us write 

Zit-T) = Zit)-TZit) + ^Z{t)... 

Substituting in Eq. (|^), the following two equations for 
the first two orders can be written. 

0(r): 

Z{t) = [(1 + iu- I Z{t) \^)Z{t) - KZ{t)] /(I - Kt), (27) 

Z{t) + bZ{t) -c{l- K + iuj- I Z{t) \^)Z{t) = 0, (28) 

where b = 2(1 - Kt)/Kt'^ and c = 2/Kt'^. It can be 
checked easily that the asymptotic periodic solutions are 
given by 



0(t) : Z{t) = 



(29) 



Oir'): Z{t) = 



1-K 



2 ll/2 



2(1 -Xt)2 



e"^'\ (30) 



where uj' = uj/{l — Kt). The effect of small r is evident 
on the frequency and the amplitude: both show a slight 
rise. 
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FIG. 4. The frequencies and amplitudes of the multiple periodic states are plotted as a function of / = tuj/{2'k) in (a) and 
(b) respectively aX K = u = 30. At large values of r, the frequency suppression is algebraic and is proportional to 1/r. In both 
the plots the unstable portions are indicated by dashed lines. 
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III. NONLINEAR FEEDBACK EFFECTS 

In this section we include the nonhnear feedback term 
{K2 ^ 0) and examine the dynamics of the hmit cycle 
oscillator in the presence of the complete feedback term 
/ as given in (2). We then have 



Z{t) = (a + iuj- I Z{t) \^)Z{t) - KiZ{t - r) 



(31) 



The transformations Z{t) e™Z{t) and K2 —K2 
leave the equation unchanged. The bifurcation diagrams, 
thus, appear symmetric about K2 — axis with the ex- 
ception that the orbits are rotated by an angle of tt about 
the origin. 



A. r = 

In order to distinguish the additional effects arising 
from the nonlinear feedback term we first turn-off the 
time delay and examine the dynamics for t = 0. These 
solutions also correspond to the special class of solutions 
with time delay when the solutions have a periodic re- 
currence with a period r. These are generally termed as 
phase trapped solutions in the literature. Let Z = X+iY. 
The evolution equations are 

X = {a-X'^ - Y'^)X -ujY - KiX - K2{X^ - Y^), (32) 
Y = ujX + [a - X"^ - Y'^)Y - KiY - 2K2XY. (33) 



1. Fixed points 

The fixed points of the system are found by solv- 
ing X — and Y ~ simultaneously. In con- 
trast to the earlier case of a linear feedback drive, 
there are now three fixed points of the system in- 
cluding the origin. In Table 1, the fixed points and 
expressions for the eigenvalues of linear perturbations 
around these fixed points are listed using the definitions 

TABLE L 



Symbol 




index 
.) 


Fixed point 

{x;,Y*) 


Eigenvalues 
A 


D 
FPi 

FP2 





1 

2 


(0,0) 
(-# + ^/G,^) 

(-^-^.^) 


—Ki ± iiu 
Ki ± Tl 
Ki ± T2 


G = 


\ 2 


? - Ki 


K, = K, 


— a and Ti — 



{X*'^ + Fj*^)^ — 07^. As can be seen, the origin ex- 
ists as a fixed point for all the values of the parame- 
ters Ki,K2,a and lo. However, it is stable only when 
Ki > 0, i.e. when Ki > a. The other two fixed points 



FPi and FP2 exist in the region where G > 0. This 
condition gives rise to K2 > 7(^1) and K2 < — 7(A'i), 



where j{Ki) = V^y Ki + y + uP-. This region over- 
laps with ki > 0. Fig. I shows the bifurcation diagram. 
The regions marked with D, FPi and FP2 are the regions 
in which the corresponding fixed points are stable. The 
region D is bounded on the right hand side by Ki ~ a. 
The region FPi is bounded by the curves Ki = and 
K2 ~ 7(^1). The region FP2 is bounded by the curves 
Ki = and K2 — —j{Ki). The movement of these fixed 
points and the corresponding eigenvalues are shown in 
Fig. |. 













Trans-critical 


PS Saddle-Node 

— ---A 


Hopf D 





FIG. 5. The bifurcation diagram of Eq. ( ^l| ) for r = in 
{Ki, K2) plane for a = 1, to = 1. D is the region of amplitude 
death, FPi and FP2 axe the regions of stability of the fixed 
points FPi and FP2 respectively. The periodic solution is 
stable in the region PS. 



As Kg is varied from Kg c to 10 
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FIG. 6. The movement of the fixed points of Eq. ( pl[ ) for 
T — {Ki = —1, a = 1, and = 1) is shown in (a) as K2 is 
varied from the critical value if 2c ~ 0.6871 to 10. The behav- 
ior of the real parts of the eigenvalues for the corresponding 
fixed points are plotted in (b) and (c) respectively and show 
that FPi is stable and FP2 is unstable in this case. 
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FIG. 7. The variation of the period of the hmit cycle in 
the region PS of Fig. ^. as a function of K2 for a fixed value 
of K\ = —1. As K2 is increased from to '){K\) the period 
increases from 2ii/uj to infinity with a scaling proportional 
to 1/ {K2c — K2) and leading to a saddle-node bifurcation. 
The inset displays on a log scale the behavior close to the 
critical region. The solid curve is a plot of 3/\/ K2c — K2 ■ 



2. Periodic orbit 



Eqs. (^2|) and ( |33D have a stable periodic solution in 
the region — 7(^1) < K2 < ^{Ki). The periodic orbit 
develops a saddle {FP2) and a node {FPi) as if 2 is in- 
creased above 'y{Ki). Similarly as K2 is decreased below 
—^{Ki) a node {FP2) and a saddle {FPi) are born. The 
periodic orbit is a circle for K2 — and the oscillator 
moves linearly on the circle with a frequency uj. As K2 
is increased or decreased the period of the limit cycle in- 
creases and tends to infinity on the curve K2 = ±"f{Ki) 
giving birth to a saddle-node point. In Fig. ^ we have 
plotted the period of the limit cycle as a function of K2 
for Ki = —l,a — l,uj — 1. The inset of this figure also 
shows the scaling of the approach to infinite time period 
which in this case is the typical inverse square root scal- 
ing of the saddle-node bifurcation phenomenon. As can 
be seen from Table 1, one of the eigenvalues of the Ja- 
cobian of Eqs. (32) and (|3^) is zero at the critical value 
of A'2. The existence of the saddle- node bifurcation can 
thus be rigorously established by reducing the dynamics 
to the center manifold Q]. 

Phase slips. As the magnitude of K2 is increased from 
to zt'^{Ki) the evolution of the phase of the oscillator 
takes on a nonlinear character. The parabolic nature of 
9{9) at its minimum and its sinusoidal nature for other 
values of 6 indicate the existence of a nonuniform mo- 
tion of the oscillator from to 2tt. The motion becomes 
jerky as the parameter is increased to bring a tangency 
at which point 9 — 0. The time period becomes singular 
at the point of tangency. 




FIG. 8. Phase slips of the limit cycle oscillator for 
T = 0,Ki = —l,a — l,u) — 1. The phase of the limit cy- 
cle oscillator evolves linearly on the circle when there is no 
nonlinear feedback. With increasing K2 one clearly sees the 
phase slips suffered by the oscillator. 

It begins to spend more and more time in the vicinity 
of 6* = 6** and less and less time at other values. Close to 
the curves ±7(^1) the oscillator shows clear phase slips 
of 27r as seen in Fig. |[ This phase slip ends with phase 
quenching on the boundary of ±j{Ki). In this state the 
phase of the oscillator is a constant and such a state is of- 
ten referred to as phase death. The phenomenon of phase 
slips plays an important role in the process of transition 
towards a synchronized state in a system of a large num- 
ber of coupled limit cycle oscillators and has been the 
subject of some recent investigations mW. 



B. T > 



We now return to our full feedback model and include 
the effect of finite time delay on the dynamical charac- 
teristics of the oscillator. Since we now have an infi- 
nite dimensional system, the phase space dynamics can 
in general be quite complicated and changes significantly 
as a function of the time delay parameter t. For our fur- 
ther analysis and description, we write down the Eq. ( |3l| ) 
in polar coordinates. 

r{t) = [a - r^{t)]r{t) - Kir{t ~ r) co%{Q{t - r) - e{t)\ 

-K2r^{t-T) co%\lQ{t - t) - Q{t)\ (34) 



-K- 



r{t) 
r^{t-T) 



r{t) 



sm[e{t - r) - e{t)] 

sin[26'(i- r) - 6l(i)]. (35) 
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FIG. 9. Bifurcation diagram of Eq. (|l|) for T = 0.5, a = 1, 
u = 1. The areas marked FPi and FP2 refer to stability 
domains of these fixed points, the region W — 1 sustains limit 
cycle orbits encircling the origin (winding number 1), the re- 
gions BSi,BS2 contain combinations of phase reversing, ra- 
dially trapped or spiraling solutions in a birythmic existence, 
and the regions marked CR- 1 ,CR- II ,CR- III ,CR- IV 
exhibit chaotic behavior interspersed with periodic windows. 

In the following we describe several different solutions 
which are found numerically (depicted in Fig. |^ and the 
following figures) and provide explanation in terms of 
the behavior of the amplitude and the phase. These 
solutions include birhythmicity, phase reversals, phase 
slips, radially trapped orbits, phase jumps, and spiral- 
ing orbits. For all the numerical integrations, unless 
stated otherwise, we use a constant history function: 
Z{t) = Z{0), t G [— T, 0) and a constant step size (h) 
of integration. Several values of h were tested for the 
accuracy of the results. Issues relating to the numerical 
integration of delay differential equations can be found in 
p2[ . In all our later results we set the parameters a — 1 
and uj — 1. 

Before we be gin our description of the solutions to 
Eqs. ( |3^ ) and (|35|), we wish to point out an impor- 
tant connection between the present model and some 
of the other models that have been studied in the 
literature. Under the assumptions that Ki = 0, 
I 9{t - r) - e{t) \« 1, and r{t) = r{t - r) = 1, Eq. (||) 
yields 

e{t) ^LU-K2 sme{t~T) 

- K2 [e{t - t) - e{t)] cos e{t - t). (36) 

The above equation with the first two terms on the right 
hand side describes a first order phase locked loop with 
time delay It also has the structure of the equa- 

tion that can be used to model certain visually guided 
movements of biological limbs |pG] . 



1. Stability of the fixed points 

We proceed first to carry out a linear perturbation 
analysis of Eq. (|3^) around the three fixed points and 
then construct the bifurcation diagram for fixed values 
of r through detailed numerical investigations. The anal- 
ysis of the stability of the origin remains unchanged by 
the presence of the K2 term since its contribution van- 
ishes in the linear limit. So its behavior can be discerned 
from the previous section (e.g. see Fig. |l|(a)), where we 
saw that the amplitude death region shrinks and moves 
to the right of the curve Ki = as r is increased from 
and vanishes at a certain critical value. Depending 
on the strength of lu, it may reappear. For the present 
study we choose uj = 1 and r = 0.5, which is much less 
than the intrinsic time period 27r. For this value the 
amplitude death region disappears and the origin is al- 
ways unstable. The fixed points FPi and FP2 are again 
the same as discussed in the previous section, but their 
stability properties are now significantly modified due to 
finite time delay effects. Let Z* — {X* + iY*) be one of 
the non-zero fixed points of the system. A linearization 
about Z* of (|l|) yields 

i(t) = {a + iuj-2 \ Z* p)z(<) 

-Z*^z{t) - {Ki + 2K2Z*)z{t-T), (37) 

where z = Z — Z* and z is the complex conjugate of 
z. Write z{t) — x{t) + iy{t) and assume each compo- 
nent to vary as e^*, where A is the eigenvalue of the lin- 
earized matrix, M2, whose eigenvalue equation is given 
by det(M2 — XI) — 0, and can be written as 

{A - Be-^^ - \f - 2(X*2 + r*^)(A - Be-^^ ~ A) 

+ {2K2Y*e-^^ -ujf ^{), (38) 

where A = a- (X*^ + Y*"^) and B = Ki + 2K2X* . The 
evolution of A in the complex plane as Ki and K2 are 
varied decides the stability of the fixed points. The re- 
lation between Ki and K2 for the critical boundary is 
not very transparent from the above equation. We use 
Eq. ( 13^ ) to determine the boundaries of the regions of 
these fixed points numerically and present the bifurca- 
tion diagram in Fig. ^. The regions labelled FPi and 
FP2 represent the regions of stability of the respective 
fixed points. When K2 is small and one is in the region 
where the origin is unstable but FPi and FP2 are stable, 
the limit cycle orbit preserves its identity and the results 
are similar to those discussed in earlier sections. As K2 is 
increased in magnitude, however, several interesting new 
orbits depicted in the bifurcation diagram appear. We 
discuss these features below. 



2. Periodic solution 

The simple periodic orbit (limit cycle) exists for all 
small values of K2 in the region marked as W-1 in Fig. ||. 
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As K2 is increased this stability can be lost through a va- 
riety of bifurcations (e.g. saddle node bifurcation, bista- 
bility or chaos). The most interesting effect of a finite de- 
lay time appears to be on the nature of the saddle-node 
bifurcation. Although the period of the limit cycle be- 
comes infinite at the bifurcation point, the scaling behav- 
ior and the nature of the orbit dynamics near this point 
are significantly different from those of the t — case. To 
illustrate this point, we have plotted in Fig. |lO|(a) a set 
of ^ vs. 9 curves for t = 0.3, i^i = as K2 is successively 
increased towards the critical value. Note that the na- 



ture of the 9 curve is not parabolic near the critical point 
and has distinct asymmetries. Further, d9/d9 develops a 
discontinuity at this point as shown in Fig. |l^(b). Since 
9 goes to near the critical point, the period of the limit 
cycle keeps increasing as K2 approaches the critical value 
as shown in Fig. ^(c). The discontinuity of d6/d9 is a 
significant deviation from the standard conditions of a 
saddle-node bifurcation |Q and brings about the change 
in the scaling behavior. To analytically estimate the time 
period close to the bifurcating point, the 9 curve, in the 
vicinity of the critical point, can be approximated as 




■ ♦ 45 io-« 10-= 10-" 

Kgc - K2 

FIG. 10. The nature of (a) d and (b) dO/dO as a function of 9 close to the bifurcation point, and (c) the corresponding 
stretching of the time period of the periodic orbit, (d) A semilog plot of the numerically obtained time periods T vs. {K2c — K2) 
(shown as circles) which indicates a logarithmic scaling of the singularity, viz. T — —3.3665 log{K2c — K2) — 1.6772 (solid line). 
Here K2C. = 0.923488, /S'l = 0, a = 1, = 1 and r = 0.3. 




FIG. 11. Periodic orbit losing stability by colliding with the saddle, (a) periodic orbit just before the critical value of K2, 
(b) corresponding evoluation of S as a function of 6. (c) shows the bifurcation diagram as K2 is increased. R is the maximum 
value of the asymptotic state of the oscillator. The parameters are t = 1.0, = —1, a = 1, and uj = 1. 
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—bix + ci, 

b2X + Ci, 



if X < 
if a; > 



(39) 



where x — ^ 9* , 9 is minimum at 9* for a given value of 
K2, and bi,b2 and ci are all positive constants, (ci = 
when K2 = Ji^2c)- The time period across a thin region 
26 around a; = is now given by 

T'. - . . , 

_s —bix + ci Jq b2X + d 



1 



1 



logci 



bi b2' 

+ log ((61,5 + cif'^^ {b2S + ci)i/^^) (40) 

As seen from the first term in the expression above, the 
time period goes to infinity as ci — > with a logarithmic 
scale. This is in contrast to the inverse square root scal- 
ing for the r = case. The above scaling agrees quite 
well with our numerical results as shown in Fig. 10|(d). 
The actual mechanism of the loss of stability of the pe- 
riodic orbit is through a collision with the saddle point. 
This is illustrated in Fig. pi] which is plotted for Ki = —I 
and T = 1.0. The curves BC and BA are the stable (FPi) 
and unstable {FP2) branches, and they are the same as 
those depicted in Fig. ^(a). As seen in the (X, Y) plane in 
Fig. ^(a) , the orbit collides with an existing saddle point 
and is kicked to the stable node. Fig. [lT|(b) shows the 
corresponding behavior of 9 vs. 9 just before the critical 
value of if 2- Notice that its intersection with the ^ = 
line is non-parabolic. There are also additional intersec- 
tions which however are not fixed points since r 7^ at 
these points. These loops have significance fo r the p hase 
reversal orbits which we discuss in Section [II B 3 , In 
Fig. pl](c), R is the maximum amplitude of the asymp- 
totic state of the orbit at each K2 ■ 



3. Phase reversals 

While still inside the region W-1, the periodic orbit 
can show a reversal of its phase with time with the same 



period as that of the orbit; the curve 9 {9) develops a 
fold as K2 is increased. Depending on the value of Ki 
there can be one or more than one folds. The acquired 
additional loop as shown in Fig. |2|(a) is not around the 
origin; the winding number of the periodic orbit con- 
tinues to be one. However the phase of the orbit as 
measured from the origin undergoes a reversal in that re- 
gion (Fig. |l^(b)), in contrast with the phase slip behavior 
discussed earlier. This is purely the effect of the time de- 
lay appearing in the nonlinear feedback term. The phase 
reversing orbits are prominent at the boundaries of the 
W-1 region. Such periodic orbits for a single oscillator 
have been observed for externally driven systems and 
basically arise due to the excitation of higher harmonics 
from the resonant interaction of the external driver with 
the basic oscillator They cannot exist for a single 

autonomous oscillator in the absence of time delay due 
to dimensional constraints. However the presence of time 
delay in our autonomous model increases the dimension- 
ality of the system and hence permits the existence of 
such orbits. This is a new kind of orbit which does not 
seem to have been noticed or discussed before in the con- 
text of large systems of coupled oscillators and it would 
be worthwhile looking for their existence in such systems. 



4- Effect of time 



on phase 



The phenomenon of phase slips continues to exist for 
non zero values of time delay as well. To study the effect 
of time delay on the phase slips, we examine Fig. |8| and 
choose the parametric values corresponding to the phase 
slips shown by the bottom most curve, i.e. for K2 — 0.68, 
and introduce finite time delay. The results are shown in 
Fig. |l^. Time delay has the effect of reducing the sharp- 
ness of the phase slips and, at the same time, increasing 
the angular speed as seen by the slope of 9. For longer 
time delays, (e.g. for r = 2.0 in the figure), phase shps 
are less evident and in fact the orbits exhibit chaos. 
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FIG. 12. Phase reversal, (a) A phase reversing orbit of the limit cycle in the X — Y space for K2 ~ 4.0. (b) Temporal 
evolution of the phase as measured from the origin as a function of K2 as K2 is increased from to 4 in steps of 0.5. (c) The 
multivaluedness of 6(6). The other constants are Ki = 1.4, a = l,tj = 1, and r = 0.5. 
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FIG. 13. The efTect of time delay on phase slips. {K2 = 0.68, Ki = —1, a = 1, and w = 1.) 




FIG. 14. Radial trapping, (a) A radially trapped orbit of the limit cycle in the X — Y space. ★ and o show two points 

separated in time by r. (b) The bifurcation diagram plotted as a function of K2- R is the maximum amplitude of the asymptotic 
state of the oscillator, (c) The bifurcation diagram in (t, ^^2) space. The fixed point FPi loses stability in a Hopf bifurcation 
as the curve is crossed to the right side. The other parameters are r = 0.5, K\ = —0.6, o = 1, and a; = 1. 
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FIG. 15. A spiraling orbit. The orbit of the limit cycle in the X — Y space (a), the temporal evolution of the individ- 
ual components X and Y (b), the amplitude r (c), and that of the phase (d) as measured from the origin are plotted for 
K2 = -3., Ki = 0.2, o = 1, w = 1, and r = 0.5. 
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5. Radially Trapped Solutions 

As we cross the boundary of FPi , the node FPi loses 
stabiHty in a Hopf bifurcation as K2 is increased result- 
ing in a periodic orbit encircling FPi. For negative K2, 
the corresponding node is FP2 . These orbits have wind- 
ing number zero around the origin and exist in the re- 
gions BSi and BS2 bordering those marked FPi and 
FP2. We call them radially trapped solutions because, 
viewed from the origin, they are restricted to a region 
of phase space and seem to oscillate within a restricted 
physical space. One such radially trapped orbit is shown 
in Fig. |l4|(a) with the corresponding bifurcation diagram 
in Fig. |lj(b). The oscillator moves clockwise. Increased 
delay requires only weaker feedback to destabilize the 
node but the value saturates as shown in Fig. |l^(c); the 
curve shown is plotted using Eq. (|3|) and is verified nu- 
merically. 

6. Spiraling solutions and phase jumps 

Another interesting periodic orbit that we find is shown 
in Fig. |l^(a) and can best be described as a spiraling or- 
bit since the amplitude of the limit cycle (in the {X, Y) 
space) first spirals out and then comes back to its origi- 
nal state. The phase changes have a step like character 



and often exhibit large jumps. Such orbits exist in a thin 
region near the radially trapped solutions. 



7. Birhythmicity 



The phase reversing and the radially trapped solu- 
tions can also coexist in the same parameter regime (e.g. 
K2 ~ —1.6, Ki — —0.2, a = 1, a; = 1, T = 0.5 as shown 
in Fig. |l^(a)) and thus exhibit a birhythmic behavior. 
In fact, birhythmicity appears to occur also for the spi- 
raling and trapped solutions and is spread out over a 
large region of the parameter space inside BSi and BS2 
in Fig. |9| Birhythmicity is a common phenomenon in 
many biological cell models, and is currently the subject 
of many studies |^^. The switching between the two 
states can take place by the slightest perturbation to the 
initial states (-'^"(O), Y{0)). A plot of the basins of attrac- 
tion for the two states of Fig. |l^(a) is shown in Fig. p^(b) 
in the space of the initial conditions {X{Q),Y{Q)). The 
dotted region is the basin of attraction for the phase re- 
versing solution and the white region is that for the ra- 
dially trapped solution. To generate this figure the initial 
conditions chosen were X{9) = X{t — 0),9 € [— t, 0), 
and y(6i) = r(i = 0), e [-r, 0). 



> 




(b) 



O 



raj'.' Jil,\ ■''a 

m- • M a . SBBBt \ "iiiiim. "-u. ""SSS^m- 




X(t=0) 



FIG. 16. Birhythmicity between two periodic orbits, (a) Radially trapped and the phase reversal solution coexisting for 
K2 = —1.6, Ki = —0.2, a = 1, tj = 1, and r = 0.5. (b) A plot of the basins of attraction for these two states in the space of the 
initial conditions {X{0),Y{0)). The dotted region is the basin of attraction for the phase reversing solution and the white region 
is that for the radially trapped solution. To generate this figure the initial conditions chosen were X(6) — X(t = 0),6 £ [— r, 0), 
and y(6») = Y{t = 0),e» G [-t,0). 
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FIG. 17. The radially trapped orbit of the limit cycle can 

undergo a period doubling sequence to the chaotic state. 
Sonic of the periodic orbits with (a) Pcriod-1, K2 ~ 1.6, (b) 
Period-2, K'2 = 2.1, (c) Period-4, K2 = 2.13, (d) Period-8, 
K2 = 2.18, (e) Period-16, K2 = 2.186, and (f) a chaotic orbit, 
K2 = 2.6 are shown. The other constants are Ki = —0.2, 
a = 1, a; = 1, and r = 0.5. The scale of (d) is twice that of 
the other plots. 




FIG. 18. The phase reversing orbit can bifurcate into a 
period-3 orbit leading to a chaotic attractor. These typi- 
cal states are shown for (a) K2 = 2.2, (b) K2 = 2.5, (c) 
K2 = 2.712. The other constants are Ki = 4.2, o = 1, a; = 1, 
and T = 0.5. 
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FIG. 19. The spiraling orbit period doubles once and then 
goes to a chaotic state as typically shown in the sequence, 
(a) K2 = 3.0, (b) K2 = 3.121, (c) K2 = 3.128. The other 



constants are Ki = 0.2, a = 1, 



1, and r = 0.5. 




FIG. 20. Typical bifurcation diagrams generated by 
Poincare section techniques showing (a) phase trapped, (b) 
phase reversing, and (c) spiraling solutions. The plots show 
the value of (a) X in the plane Y ^ 0, Y < 0, Ki = -0.25, 
(b) Y in the plane X = -lA, X < 0, Ki = 4.3 and (c) X in 
the plane Y = -2., Y < 0, Ki = 0.2. The other constants 
are a = 1, w = 1, and t = 0.5. 
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FIG. 21. Bifurcation diagram as a function of r. R is the maximum value of the amplitude of the oscillator as a function of 
r for K2 = -3.0, and = 0.2, a = 1, and = 1. The initial conditions used are X(0) = 0.123, F(0) = 0.456. The orbits at 
various points that are numbered are shown in the inset. The orbit at 11 and 13 is the spiraling orbit with period-2. 



8. Routes to Chaos 



The nonlinear model equation (|31|) also exhibits re- 
gions of temporal chaos and we have investigated this 
phenomenon in some detail, particularly with regard to 
routes to chaos. The route to chaos appears to be a 
strong function of the parameters Ki , K2 and r and the 
type of periodic orbit existing in a particular parameter 
domain. There is evidence of several distinct routes to 
chaos in the present model system. For example, the ra- 
dially trapped solutions appear to follow the period dou- 
bling route to chaos as shown in the sequence of plots 
in Fig. The phase reversing solutions appear to go 
to the chaotic state in several different ways. The limit 
cycle can make a transition from a simple orbit to a quasi- 
periodic state. It can also period double. Or a period I 
state can make a transition to a period-3 orbit. Such a 
transition is shown in Fig. We find a large window 
of Pcriod-3 orbits starting from K2 — 2.283 for the pa- 
rameters given in Fig. |l8|. The spiraling orbits appear 
to undergo a single period doubling bifurcation and then 
a sudden transition to chaos. An example is shown in 
Fig. |l^. Finally, in Fig. ^ we display the detailed bifur- 
cation diagrams (obtained from Poincare section plots) 
corresponding to these scenarios. We would like to add 
a word here about the numerical precautions followed 
in generating these diagrams. Since many of these at- 
tractors can coexist in a birhythmic state, it is necessary 



to adopt high accuracy and care in tracking any one of 
them. Likewise one also needs to get beyond some long 
transients that the system can exhibit for some particular 
initial conditions. Using these precautions, in Fig. |2l|, we 
have given a detailed bifurcation diagram as a function of 
the time delay parameter. In this, the maximum ampli- 
tude of the oscillator is plotted for Ki = 0.2, K2 — —3.0 
at various values of r. The various numbers on the curves 
stand for the kind of orbit that is found in that region (a 
phase portrait of the orbit is also depicted for reference). 
The initial conditions chosen have been X{Q) — 0.123 
and r(0) = 0.456. 



IV. SUMMARY AND DISCUSSION 



We have studied the dynamics of a single Hopf bi- 
furcation oscillator (the Stuart-Landau equation) in the 
presence of an autonomous time delayed feedback. The 
feedback term has both a linear component and a sim- 
ple quadratic nonlinear term. Using a combination of 
analytical methods and numerical analysis, we have in- 
vestigated the temporal dynamics of this system in var- 
ious regimes characterized by the natural parameters of 
the oscillator (e.g. its frequency cu, linear growth rate 
a), strengths of the feedback components {Ki, K2) and 
the time delay parameter, t. Our principle results are 
presented in the form of bifurcation diagrams in these 
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parameter spaces. These reveal a rich variety of tempo- 
ral behavior including time delay induced stabilization of 
the origin, multiple frequency states, frequency suppres- 
sion, phase slips, saddle node bifurcations, and chaotic 
behavior. In addition, some of the periodic orbits ex- 
hibit novel behavior such as birhythmicity, phase rever- 
sals, radial trapping, spiraling oscillations in amplitude 
space. Some of these can be understood by the stabil- 
ity of the fixed points or the loss of single valuedness of 
amplitude evolution equation. 

One of the attractive features of these results is that 
many of them have been observed in the collective be- 
havior of larger systems such as the Kuramoto model or 
the amplitude versions of the Kuramoto model. This has 
been one of our major motivations for constructing this 
model - as a sort of paradigm to obtain and investigate 
these states in a simple manner. The feedback terms not 
only model the collective drive that a single oscillator 
feels in a larger system, but also incorporate time delay 
in an autonomous manner. Time delay frees the dimen- 
sional constraints of the system (the system is essentially 
oo -dimensional) and this might be the reason why its 
temporal dynamics resembles so much that of larger di- 
mensional systems. Our results may therefore be useful 
for gaining better insights into the behavior of such large 
systems. As an example, the large window of — 3 or- 
bits of our model appears to have a correspondence to 
the transition region between the incoherent and chaotic 
states of coupled limit cycle oscillators. We have found 
that the oscillators in the wings of the frequency dis- 
tribution of a large collection of oscillators in the mean 
field model begin to acquire the 14^ — 3 temporal states in 
that region and play an important role in the transition 
mechanism |Q. A detailed understanding of their dy- 
namical behavior can thus help us in addressing some of 
the outstanding problems in this area, such as the nature 
of transition between low dimensional chaos and turbu- 
lence. Our model can also find more direct applications 
in simulation studies for feedback control of individual 
physical, chemical or biological entities that have the ba- 
sic nonlinear characteristics of our Hopf oscillator, such 
as in single mode semiconductor lasers, relativistic mag- 
netrons, chemical oscillations, and biological rhythms in 
single nerve cells. In fact, the basic Stuart-Landau equa- 
tion is mathematically related to the well known van der 
Pol oscillator equation from which it can be derived by 
a suitable time averaging. Our oscillator model with the 
quadratic nonlinearity can likewise be derived from a van 
der Pol type equation which in addition has a nonlinear 
Mathicu like term, i.e., it is an equation of the form 

X + n {x'^ - l)x + {1 + ex sin(t)) x = 0. (41) 

Such a nonlinear equation can physically represent para- 
metric excitation of relaxation oscillations and can be 
used to model a number of physical or biological systems. 
It may be possible in such simple systems to then seek 
experimental verification of some of the novel temporal 
states displayed by our model. 
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